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Abstract. Extremal Optimization (EO), a new local search heuristic, is used to approximate ground states 
of the mean-field spin glass model introduced by Sherrington and Kirkpatrick. The implementation extends 
the applicability of EO to systems with highly connected variables. Approximate ground states of sufficient 
accuracy and with statistical significance are obtained for systems with more than A'^ — f 000 variables using 
it J bonds. The data reproduces the well-known Parisi solution for the average ground state energy of the 
model to about 0.01%, providing a high degree of confidence in the heuristic. The results support to less 
than 1% accuracy rational values of lo — 2/3 for the finite-size correction exponent, and of p = 3/4 for the 
fluctuation exponent of the ground state energies, neither one of which has been obtained analytically yet. 
The probability density function for ground state energies is highly skewed and identical within numerical 
error to the one found for Gaussian bonds. But comparison with infinite-range models of finite connectivity 
shows that the skewness is connectivity-dependent. 
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1 Introduction 

The Sherrington-Kirkpatrick (SK) model has provided 
a rare analytic glimpse into the nature of frustrated spin 
glasses below the glass transition. It extends the notion 
of a spin glass on a finite-dimensional lattice introduced 
by Edwards and Anderson (EA) |2| to infinite dimen- 
sions, where all spin variables are infinitely connected and 
mean-field behavior emerges. In this limit, analytically in- 
tractable geometric properties of the lattice submerge. 
Consequently, the SK model simply establishes mutual 
bonds between all variables. Many features of this highly 
connected model have become analytically accessible with 
Parisi's replica symmetry breaking (RSB) scheme j2i- Only 
recently have RSB models with long-range but finite con- 
nectivity been analyzed successfully An comparable 
treatment of EA is still missing. 

The SK model remains a topic of current research 
EIIZI- For one, its mathematical challenges, leaving cer- 
tain scaling exponents as-of-now intractable, continue to 
inspire new theoretical approaches [H|. Furthermore, as 
scaling arguments |51llU| for EA suggest an entirely dif- 
ferent picture, the fundamental question to the relevance 
of mean-field theory for any description of realistic sys- 
tems at low temperature remains unanswered. 

The challenge of the SK model is exemplified by the 
fact that it is an NP-hard problem to find the ground 
state of its instances |31 . Unlike in a spin model of ferro- 
magnetism, in which couplings Jij — 1 always try to align 
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neighboring spins, in a spin glass model like SK or EA, 
each spin is frustrated by a competition between randomly 
drawn, aligning and anti-aligning couplings (say, Ji,j = 
±1) to its neighbors. As a result, its potential energy land- 
scape is characterized by a hierarchy of valleys within val- 
leys with a number of local minima growing expo- 
nentially in the system size . Since its low-energy land- 
scape features prominently in its low-temperature proper- 
ties, even numerical insights have been hard to come by. 
Some earlier work in this area has been focused on gradient 
descent |[12,13 or Simulated Annealing algorithms |14| . 
extrapolations to low temperatures from perturbative ex- 
pansions near the glass transition |15) . or on exact meth- 
ods to enumerate low- lying energy values ^H]. And even 
with the most sophisticated methods, like genetic algo- 
rithms (GA), accurate approximations have been limited 
to system size of < 300 PMT) . 

Here, we propose an alternative optimization proce- 
dure, based on the Extremal Optimization (EO) heuris- 
tic |17lll8j . Our implementation of EO 10.^ is extremely 
simple and very effective, allowing to sample systems of 
sizes up to N 1000 with sufficient accuracy and statis- 
tical significance. This approach produces results that not 
only verify previous studies by independent means, but 
also improve the accuracy. Previous studies [31] suggest 
that the fluctuation exponent of the ground state energies 
p is near to 3/4, excluding an earlier conjecture of 5/6 j^U] 
1^ . Here, we double the size of the scaling regime to find 
p = 0.7500(29). These results strongly support analytical 
arguments by Refs. |22llK] in favor of p = 3/4, assuming 
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that such an exponent in a solvable model should be a 
simple rational number. 



2 EO Algorithm 

Our implementation of t-EO proceeds as follows |17II18| : 
Assign to each spin variable Xi{= ±1) a "fitness" 



N 



Aj — Xi ^ ^ Ji^jXj^ 

i. e. the (negative) local energy of each spin, so that 



(1) 



N 



H 



2VN 



(2) 



is the familiar Hamiltonian of the SK model. For general 
bond matrices Jij , such as those drawing from a continu- 
ous Gaussian bond distribution with varying bond- weights 
attributed to different spins, more refined definitions of 
should be used 17, 25 . Here it is conceptually and compu- 
tationally most convenient to draw discrete bonds J from 
{ — 1,+1} with equal probability, such that (J) = and 
(J2) = 1. 

A local search with EO TT^ ideally requires the ranking 
of the fitnesses A^ from worst to best before each update. 



Ai7i < A 



< 



< A 



(3) 



where i — Ilk indicates spin Xi as having the fc-th ranked 
fitness. At each update, one spin of low fitness is forced 
to change unconditionally. Since EO does not converge to 
a specific configuration, it outputs the best-found after a 
certain number of updates. 

Following Ref. jl7j . it is most expedient to approxi- 
mately order the Xi in Eq. ^ instead on a binary tree of 
depth 0(log2 N) with the least-fit spins ranking near the 
root. Unlike for sparse bond-matrices 18 , flipping one 
spin also changes the fitness of all other spins, albeit by 
a small amount, AXi/Xi = 0{1/N). To avoid the cost 
of O(logA^) for re-ordering the entire tree each update, 
a dynamic ordering scheme is used here: All A^ are re- 
evaluated, but the tree is parsed only once, node-by-node, 
starting at the root. The fitness on the current node is only 
compared with its two sub-nodes and exchanged, iff its fit- 
ness is better. In this way, a newly improved fitness can be 
moved away from the root several times, but newly worse 
fitnesses move at most one step towards the root. Yet, a 
spin which suddenly attained a low fitness would move 
to the root at most within 0(log2 N) updates. Hence, re- 
ordering of fitnesses occurs faster than mis-orderings can 
escalate because AX/X ^ 1. 

In a r-EO update, a spin is selected according to a 
scale-free probability distribution P{k) ~ k^"^ over the 
ranks k £ {1, . . . , N} in Eq. ||2Jl. Since the ranking here 
is not linear as in Eq. but on a tree, a level ^, < 
I < [log2(7i)J is selected with probability ~ 2^^^^^^', and 
one randomly chosen spin on the Z-th level of the tree 
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Fig. 1. Plot of the average best energy per spin found by 
EO as a function of the parameter r. For each system size 
N, a set of test instances were created and optimized with 
r-EO, each for update steps. Each data point represents 
the average over the best-found energies obtained with that r. 
In accordance with Ref. 1261 . the optimal choice for r within 
the given runtime moves closer to unity slowly with increasing 
system size. Within the range of N used here, a fixed t « 1.2 
appears to be effective. 



is updated jl7| . In this manner of ranking and selecting 
from a binary tree, an ideal selection according to P{k) 
is approximated while saving 0{N) in the computational 
cost. Tests show, in fact, that the r-dependence for op- 
timal performance of this algorithms follows the generic 
behavior described in Ref. see Fig.^ EO at t — 1.2 
finds consistently accurate energies using 0{N^) update 
steps in each run, at least for N < 1000, verified by the 
fact that our data reproduces the exactly known energy 
of the SK to about 0.01%, see Fig. |21 Including the linear 
cost of recalculating fitnesses and dynamic ordering, the 
algorithmic cost is 0{N'^). Runs take between « Is for 
iV = 63 to « 20h for N = 1023 on a 2GHz Athlon CPU. 

It is not at all obvious that EO would be successful in 
an environment where variables are highly connected. So 
far, EO has only obtained good results for systems where 
each variable is connected only to 0(1) other variables for 
N —> oo. The update of a single variable hence impacts the 
extensive energy of the system only to sub-leading order, 
and only 0(1) variables need to rearrange their fitness. 
Applications of EO to highly connected systems, where 
each degree of freedom is coupled to most others over long- 
range interactions, proved unsatisfactory: For instance, in 
a continuum polymer model |27| with torsion angles be- 
tween chain elements as variables, even a minute rotation 
leads to macroscopic changes in the total energy, and al- 
most all moves are equally detrimental. In that case, cri- 
teria for move rejection are necessary, which are decidedly 
absent from EO so far. But for the SK in a update near Eq 
we estimate AE/E = E»^A</E,A» < t.^i^>'^/ >'^) ^ 
1/\/N, assuming a sum over terms with random signs. In 
fact, the ability to sustain roughly ^/N perturbations to 
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Table 1. List of our computational results to approximate 
ground state energies eo of the SK model. For each system size 
A'^, we have averaged the energies over / instances and printed 
the rescaled energies (eo), followed by the deviation a{eo) in 
Eq. Given errors are exclusively statistical. 



iV 


I 


(eo) 


cr(eo) 


15 


380 100 


-0.64445(9) 


0.0669(3) 


31 


380 100 


-0.69122(8) 


0.0405(2) 


49 


500 000 


-0.71051(6) 


0.0293(1) 


63 


389 100 


-0.71868(3) 


0.0246(1) 


99 


500 000 


-0.73039(3) 


0.01763(7) 


127 


380 407 


-0.73533(2) 


0.01468(7) 


199 


351317 


-0.74268(2) 


0.01043(5) 


255 


218 473 


-0.74585(2) 


0.00862(5) 


399 


15 624 


-0.75029(5) 


0.0061(1) 


511 


25 762 


-0.75235(3) 


0.0051(1) 


799 


725 


-0.7551(1) 


0.0037(4) 


1023 


244 


-0.7563(2) 


0.0029(6) 



the system before altering the macroscopic state may be 
one of the advantages of EO. 



3 Numerical Results 
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Fig. 2. Extrapolation of the data for (eo) in Tab.^for N oo. 
The exact result of —0.76321(3) (x) is reproduced within 
0.01% accuracy. The near-linear behavior of the fit yields a 
scaling-correction exponent of u = 2/3 to about 1%. The inset 
shows the same data, subtracted by —0.76324 and rescaled by 
N'^^^ , which now extrapolates to the amplitude of the scaling 
corrections at « 0.70(1). Despite "peeling off" layers of the 
asymptotic behavior, the data remains quite coherent, attest- 
ing to the accuracy of the EO heuristic. 



Extensive computations to determine ground state en- 
ergies per spin, eg, of about / = 5 x 10^ instances for 
N < 100 to just / « 250 instances for N = 1023 have 
yielded the results listed in Tab. ^ Note that all values 
chosen for N are odd. Using = 2* — 1 was convenient to 
ensure a complete filling of all levels on the tree ranking 
the fitnesses in Sec.|21 Subsequently, we added data at in- 
termediate values of N. For smaller N there was a minute 
but noticeable deviation in the behavior of (eo) between 
even and odd values of N, with even values leading to 
consistently lower (eo). Either set of data extrapolates to 
the same thermodynamic limit, with the same corrections- 
to-scaling exponent, but appears to differ in the ampli- 
tude of the scaling corrections. This behavior is consistent 
with the findings for even and odd-connectivity Bethe lat- 
tices |2Hi ■ (Note that even TV here implies odd connectivity 
for each spin in the SK model, and vice versa.) 

We have plotted (eg) vs. l/A^^/s Fig. El The data 
points extrapolate to —0.76324(5), very close to the best 
known Parisi energy of —0.76321(3) jJS]. All data shown 
in Fig.|21fits to the asymptotic form (eo)Ar = (eo)oo+a/A^'^ 
with a goodness-of-fit Q w 0.7. The fit gives for the ex- 
ponent for scaling corrections uj — 0.672(5), or 2/3 within 
1%. This is consistent with analytical results for scaling 
corrections obtained near Tg [211 and with numerical stud- 
ies of ground state energies 013 for the SK model, but also 
with EO simulations of spin glasses on finite-connectivity 
Bethe lattices and ordinary random graphs PU) . 

The large number of instances for which estimates of 
eo have been obtained allow a closer look at their distri- 
bution. The extreme statistics of the ground states has 
been pointed out in Ref. ^Tj and studied numerically in 
Refs. IHltZ]. Being an extreme element of the energy spec- 
trum, the distribution of eo is not normal but follows a 



highly skewed "extremal statistics" [22 • If the energies 
within that spectrum are uncorrelated, it can be shown 
that the distribution for eo is among one of only a few uni- 
versal functions. For instance, if the sum for H in Eq. (0) 
were over a large number of uncorrelated random variables 
Ai, H would be Gaussian distributed. In such a spectrum, 
the probability of finding H —^ —oo decays faster than 
any power, and ground states eg would be distributed ac- 
cording to a Gumbel distribution, |31irf] 

( X - u \x - u]\ 

gm{x)=wexp<m mexp > (4) 

I ^ [ t; J J 

with m = 1, where m refers to the m-th lowest extreme 
value. 

Clearly, in a spin glass the local energies are not un- 
correlated variables, see Eq. lO, and deviations from the 
universal behavior may be expected. In particular, these 
deviations should become strongest when all spin variables 
are directly interconnected such as in the SK model, but 
may be less so for sparse graphs. Indeed, in the SK model 
with Gaussian bonds Refs. find numerically highly 

skewed distributions for eo which do not fit to the Gum- 
bel distribution in Eq. Q for m = 1. In Fig. we plot 
the rescaled distribution of ground state energies obtained 
here for ± J bonds. The result resembles those of Ref. [7] 
to a surprising degree. In fact, a naive fit of Eq. (QJ for 
variable m to the SK-data, as suggested by Ref. [J] , yields 
virtually identical results, with m « 5. This may indicate 
a high degree of universality with respect to the choice of 
bond distribution in the SK model, or a new universality 
class of extreme- value statistics for correlated variables. In 
Fig.|31we have also included data for fc-f 1-connected Bethe 
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Fig. 3. Plot of the rescaled probability distribution of ground 
state energies using ±J bonds. Shown are the data for the SK 
model and for Bethe lattices of connectivity A; + 1 = 3 and 25 
from Ref. ^'28K The data for increasing k seems to evolve away 
from a Gaussian (solid line) towards the SK data (fc — oo), the 
latter fitted by Eq. The values obtained in the fit (dashed 
line) are u — 0.26, v = 2.23, w — 90, and m — 5.4. 
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Fig. 4. Plot of the standard deviation in the distribution of 
ground state energies eo vs the system size A'^. Asymptotic 
scaling sets in for A'' > 63, clearly favoring N~^^'^. A fit (full 
line) of these data points extrapolates to p = 0.7500(29). The 
inset shows the same data reduced by the predicted asymp- 
totic scaling, a{eo)/N~^^'^ , as a function of 1/A''^'''*. Any devi- 
ation from A'^~^'^^-scaIing would appear as divergent behavior 
for oo. Instead, the scaling corrections are well-captured, 

say, by a simple parabola in 



lattices from Rcf. 28 for fc + 1 = 3 and 25, which seem to 
suggest a smooth interpolation in k between a normal dis- 
tribution and the SK result. Hence, while the distribution 
of eo seems to be universal with respect to bond distribu- 
tion, its connectivity-dependence appears to disfavor the 
existence of a (unique) universal extreme-value statistic 
for correlated energies. 

We now consider the scaling of the standard deviations 
in the distribution of eg with respect to system size, 

'^(eo) = ^{el) - (eo)' - N~^. (5) 

where p is the fluctuation exponent. Similarly, the fluc- 
tuations of eo appear to be narrower than normal, with 
p > 1/2 in Eq. Early theoretical work ^20,21^ sug- 
gested a value of p = 5/6. More recent numerical work 
|7| instead is pointing to a lower value. Ref. |^ have ad- 
vanced an alternative argument in favor of p = 3/4, based 
on corrections in the zero-mode of the propagator due to 
fluctuations. 

In Fig. 0] the numerical results for the standard devi- 
ations in the distribution of ground state energies eo is 
shown. The asymptotic scaling for > 63 is certainly 
very close to p = 3/4. The crossover toward asymptotic 
behavior is similar to the results found for Gaussian bonds 
using a GA (see Fig. 1 in Ref. 7 ), except that the EO 
data reaches about half a decade further into the asymp- 
totic regime. A fit, weighted by the statistical error, to the 
data points in the scahng regime yields p = 0.7500(29), 
or 3/4 within 0.4%, with a goodness-of-fit Q = 1. As the 
inset of Fig.^jshows, any apparent trend towards a higher 
value [3| then p = 3/4 is easily explained in terms of scal- 
ing corrections, for instance, in powers of 1/A^^/*. 



4 Conclusions 

We have shown that the extremal optimization heuris- 
tic can be extended successfully to highly connected sys- 
tems. Results for the ground states of the SK model are 
consistent with previous studies while reaching assuringly 
larger systems sizes. These results provide more confidence 
into conjectures about as-of-yet unobtainable scaling ex- 
ponents. Comparison with data for k-\- 1-connected mean- 
field spin glasses on Bethe lattices suggest a smooth inter- 
polation in k for the extreme- value statistic of the ground- 
state energy between a Gaussian distribution for small k 
and a highly skewed Gumbel distribution with m « 5 for 
the SK model {k — + oo). 
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